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^ . Abstract 



We consider the frequency response problem and derive a posteriori 
error estimates for the discrete error in a reduced finite element model 
obtained using the component mode synthesis (CMS) method. We pro- 
vide estimates in a linear quantity of interest and the energy norm. The 
' estimates reflect to what degree each CMS subspace influence the over- 

y-^ , all error in the reduced solution. This enables automatic error control 

■ through adaptive algorithms that determine suitable dimensions of each 

(**] ' subspace. We illustrate the theoretical results by including several nu- 

merical examples. 



Keywords: Component mode synthesis; model reduction; reduced order modeling; a 
posteriori error estimation; frequency response problem. 



^ ; 1 Introduction 

t^j- ' Due to the large scale of finite element models of complex structures, it may 

l/^ . be necessary to use reduced finite element models with much fewer degrees 

' of freedom when performing frequency response analysis of a structure over 

a large range of frequencies. Having control over the reduction error in the 
O^l ' approximation is then highly important. In this paper we derive a posteriori 

error estimates for the discrete error in the reduced solution to the frequency 
response problem obtained using component mode synthesis (CMS) [6, 71 1111 HI 

El- 

The results in this paper complement the a posteriori analysis developed in 
. [13] where CMS was applied to an elliptic model problem, and the a posteriori 

analysis in |14) . where the elliptic eigenvalue problem was considered. Similar 
techniques are used and results obtained here as in the previous two publications. 
The frequency response problem does however require an explicit treatment due 
to the indefinite nature of the problem. We further present a new adaptive 
strategy suitable for frequency sweep analysis. 

Other work on error analysis for CMS class methods include results by 
Bourquin who considered the elliptic eigenvalue problem and derived a priori 
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bounds for the error in eigenvalues and eigenmodes [H0; and results by Yang, 
Gao, Bai, Li, Lee, Husbands, and Ng Q15] for the automated multilevel substruc- 
turing method [3J, who derived a criterion for mode truncation; together with 
results by Elssel and Voss [8] who showed that the same criterion guarantees 
control of the error in the smallest eigenvalue in the reduced problem. 

Previous work on frequency response analysis based on CMS include that 
by Bennighof and Kaplan [2], who developed an iterative method in which 
the response is split into two components, one component near resonance and 
one component representing the remainder of the response. The near resonant 
component is captured using approximate global eigenmodes, and the remainder 
of the response using substructure modes and iteration. A similar method was 
also proposed by Ko and Bai [T5]. Error estimates for these methods have, to 
the author's knowledge, not yet been developed. 

Research on duality based a posteriori error estimation and adaptive refine- 
ment strategies in finite element modeling has been ongoing since the 1990's. 
For a general introduction to the subject in context of finite element analysis 
we point the reader to [H [TU1 IS], and the references therein. We also refer 
to [HJ HH [T71 HH] for results that we feel are especially relevant in context of 
structural mechanics. 

The remainder of this paper is organized as follows. In Section 2 we provide 
some prerequisite material and present the frequency response problem in linear 
elasticity; in Section 3 we give an account of the Craig-Bampton CMS in a 
variational setting; in Section 4 we derive an a posteriori error estimates for the 
error in the displacements in the reduced model measured in the energy norm; 
in Section 5 we demonstrate our results in several numerical examples; and in 
Section 6 we summarize our findings. 

2 The Frequency Response Problem 

Let O be a bounded domain in R d , d = 2, 3, with boundary <9f2 = Fd U Fat, 
where n Tjy — 0, let T be a subdivision of £1 into, for instance, triangles 
(d = 2) or tetrahedra (d = 3), and let V h be the space of continuous, piecewise 
pth order vector polynomials on T defined by V = {v £ [i? 1 (f2)] d : v\-p D = 
0, v\ T e [P p (T)} d , VT e T}, where V P (T) is the space of pth order polynomials 
on element T. Let further a(-, •) be the bounded, coercive bilinear form on 
yh x yh c; ennec i by a (v,w) = 2(fie(v) : e(w)) + (/tV • v, V • w), where e(v) : 
e(w) = Yli,j=i £ ij( v ) £ ij( w )i l et (') ') denote the L 2 inner product on V h x V h , 
and let £>(•) be the bounded linear form on V h given by b(v) = (/, v) + (giy, v)p N , 
where / is a body force, and gjv is a traction force. 

The finite element frequency response problem reads: find U £ V h such that 

a(U,v) + iuj(VU,v)-u 2 (U,v) = b(v), VveV h . (2.1) 

Given a basis in V , the following matrix form of (|2.1j) is obtained: 

KU + twDU-w 2 MU = b, (2.2) 

where K is the stiffness matrix, M is the mass matrix, D is a damping matrix, 
assumed to be on the form D = aK + /3M, a > 0, /3 > 0, i.e. Rayleigh damping, 
and b is the load vector. The vector of coefficients of U is denoted by U. 
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3 Component Mode Synthesis 

Let S = be a partition of Q into n connected subdomains flj, such that 

each 0^ = \JjceK(K, for some subset /Cj c /C. Let the interface between the 
subdomains be denoted by T. An a-orthogonal decomposition 

n 

V h = @Vf, (3.1) 

i=l 

of V h associated with S and L may be constructed by letting V/ 1 = {v G V h : 
u|n\n< = 0}, i = 1, . . . ,n, and by letting 

Vj l = {£u eV h :v eV h \r}, (3.2) 

where V |r denotes the trace space of associated with T, and ft' G W l 
denotes the energy minimizing extension of a function v 6 V^lr to f2. That is, 
£i/ is defined by the problem: find £v> G y' 1 , such that 

a(Su,v) = 0, VveVj 1 , i = l,...,n, (3.3) 
Sv\ r = v. (3.4) 

A basis in each subspace Vj 1 , i = 0, . . . , n, assumed to be of dimension fcj, 
is obtained from the discrete eigenvalue problems: find (Aj, Zi) € Rx Vj 1 for 
£ = 0, . . . , n, such that 

a(Z h v) =Ai(Zi,v), VveVf, i = 0,...,n. (3.5) 

A reduced subspace V h,m C V h , where m = (m,-)^_ is a multi-index, may be 
defined by letting 

n 

where 

V? D = spanjZ^}^, i = 0, . . . , n. (3.7) 

3.1 The Reduced Problem 

Introducing the subspace V h ' m in the model we get the following reduced prob- 
lem: find U m € V h m such that 

a{U m ,v) + mj(VU m ,v) -u 2 (U m ,v) = b{v), Vt> G yh,m ^ 

Collecting the coefficients of the reduced basis functions columnwise in the ma- 
trix V" 1 , the matrix form of (|3.8p reads 

K m U™ + twD^U™ -w 2 M m U m = b m , (3.9) 

where 





= (V m ) T KV m , 


(3.10) 


D T " 


= (V m ) T DV m , 


(3.11) 


M m 


= (V m ) T MV m , 


(3.12) 


b m 


= (V m ) T b. 


(3.13) 
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Next, we turn to error estimation and derive a posteriori error estimates for 
the discrete error E = U — U m in the reduced problem. 

4 A Posteriori Error Analysis 
4.1 Preliminaries 

We begin by remarking that for the error E holds the Galcrkin orthogonality 
property 

a(E,v) + LUj(E,v)-uj 2 {E,v) = 0, Vv G v h ' m , (4.1) 

obtained by subtracting (|3.8p from the finite element formulation (|2.1[) . 

We further state the following notation, first introduced in [13] ■ Here and 
below IZi : V h — > V/ 1 , i = 0, . . . , n, denote Ritz projectors, that is, the projector 
defined by 

a(w - TZ t w, v) =0, V« e V./\ (4.2) 
and 7Z : V h —> V , denotes decomposition, such that 

n 

Kw^^KiW. (4.3) 

The operators : V- 1 — > y^. m i j i _ o, . . . , n, denote series expansion in 

V i ' mi , such that 

mi 

V? w ^V.Z,.,;Z,. r (4.4) 
i=i 

and the operator P" 1 : y' 1 — > V ,m , similarly denotes expansion in the space 
V^' m , such that 

n 

P m w = ^Pf i ^ i w. (4.5) 

i=Q 

We further let R\(w) e , i = 0, . . . , n, denote the discrete subspace residuals 
defined for w e V h by 

), V«e Vf. (4.6) 
We note here that the discrete residual R h (w) £ V h is defined by 

(R (w), v) = b(v) — a(w, v) — loj(Dw, v) + ui 2 (w, v ), \/v£V h , (4.7) 
and that the l? projection ViR h (w) of R h (w) onto V/ 1 is given by 

(R h (w)-PiR h (w),v)= WeVf. (4.8) 
Hence, R£(w) = V l R h {w) 1 and 

(R^(w),v) = (R h (w),v), VveVf. (4.9) 
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Finally, for a > we define the operators Cf : V/ 1 — > V i , £ = 0, . . . ,n, by 

(£ j a tt,«) = X;A? j (ii 1 Zy)(Zy,e) l Vti.reVJ*, i = 0,...,n, (4.10) 

3=1 

where the G , i = 0,...,n, are given by (|3.5p . We summarize some 
properties of Cf in the following lemma whose proof is straight forward, and 
may be found in e.g |13) . 

Lemma 1. For Cf , i = 0, . . . ,n, as defined in ([4.10[) . £/ie following properties 
hold: 

(CiU,v) = a(u,v), Vu,veV- 1 , i = 0,...,n, (4.11) 

||£ 4 1/2 tt|| 2 = |||«||| 2 , V«GV5 h , » = 0,...,n, (4.12) 

||(/-PrH<-^— ||£? W ||, V«eV/\ i = 0,...,n. (4.13) 

A !,m, + 1 

4.2 A Posteriori Estimate in a Quantity of Interest 

First we show a straight forward a posteriori error estimate for the discrete 
error in a quantity of interest defined by a linear functional. The quantity of 
interest may for instance be the mean stress on a part of the boundary, or the 
displacements near a point of interest. 

Let therefore H(-) — ip G V , be a linear functional on V , and let 

the goal of solving (|3.8p be to accurately approximate H(U). We introduce the 
dual problem: find <fr G V h , such that 

a(v,<!>) + LU}(v,V<f>)-uj 2 (v,<&) =H(v), G V h . (4.14) 

Choosing v = E in (I4.14p . the error H(E) may then be written 

H{E) = a(E,lZ&) + luj(VE,1Z<&) — u 2 (E,lZ$>) (4.15) 

= £(^(£7™),^*) (4.16) 

n 

= Y,{R h (U m ),K&), (4.17) 

8=0 

using (|4.9p . and using the triangle inequality, the estimate 

n 

lif^l^^KH^t/" 1 ),^*)], (4.18) 

i=0 

immediately follows. 

Now, using that {Zi,j}jLi is an orthonormal basis in V/ 1 , we have for the 
residual Rf(U m ) and the Ritz projection lZi<&, respectively 

ki 

R h {u m )= ST ■Jtdl-.Z^-Z,.,. (4.19) 

j=rrii +1 

AVt> ^ f " :<l ' N Z '' ; Z,, ; . (4.20) 



3=1 
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Thus, 



\H( E )\<± £ i^^y"),^, (421) 

j— j— mi+1 



Remark 1 . The accuracy in the reduced model of course depends on the dimen- 
sions rrii of the subspaces yh mi t % = 0, . . . ,n. Looking at the estimate (|4.18|l . 
each term 

i7 J|i = |(il h (E/ m ) 1 ^*)| 1 i = 0,...,n, (4.22) 

accounts for the contribution to the error H(E) caused by truncating the basis 
in y"/ 1 ' 1 ™* ; { = 0, . . . , n. The r]j t i may then be used as a decision basis in an 
adaptive algorithm that automatically refines the subspaces contributing the 
most to the error H(E). 

Remark 2. To obtain r)j t i, we need the residual R h (U m ) and the dual solution 
The coefficient vector R of the residual R h (U m ) is given by the equation 

MR = b - KtJ m - (DtJ m + w 2 MU m , (4.23) 

and the dual problem on matrix form reads 

K* + iD* - w 2 M* = H, (4.24) 

The coefficient vector <&; of lZi& is given by the equation 

Vf KV,*, = Vf K*. (4.25) 

where Vj is a matrix containing the coefficients of a basis in V/ 1 in its columns. 
In the case of the modal basis {Zi^} k ^ =1 , we then have V^KV^ = Aj, where 
Ai is diagonal, and we obtain 

r\ u = R T MVJ, (4.26) 

= R T MViAr 1 vfK* (4.27) 

= (b - KU m - tDU m + a; 2 MU m ) T V l Ar 1 Vf K* (4.28) 

= (b - KU m - tDU m + a; 2 MU m ) T W 4 Ar 1 Wf K$, (4.29) 

where we have used that R h {U m ) is orthogonal to V i h ' mi in the last equality, 
and introduced W,, which we assume contains the coefficients of the fcj — m, 
eigenmodes {Zij} k ' =m . +1 , together with the diagonal matrix Aj containing the 
corresponding eigenvalues. 

In practice the dual problem is approximately solved, for instance using 
a slightly larger reduced space V h ' d , where d = (rfj)"_ , and m, < d,: < fcj, 
i — 0, . . . , n, compared to what is used in the primal problem. Similarly, in the 
Ritz projections of the dual solution onto the subspaces, approximations may 
be used. Due to orthogonality it is then sufficient to project onto the spaces 
yh,di \ v h ' mi , i = 0,...,n. 
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4.3 An Energy Norm Estimate 



The following a posteriori error estimate in the energy norm ||| • ||| = y/ a(-, •) 
holds. 

Theorem 1. Let U and U m satisfy (|2.1[) and (|3.8p . respectively. Then the 
following a posteriori estimate holds for the energy norm of the discrete error 
E = U — U m in the approximation: 



\\\E\\\ < [^h + S{uj)^2T^ (4.30) 
where I\ = Ii(U m ) and I2 = l2(U m ) respectively, are defined by 

n 

h(U m ) = Y,T W R i(U m )f, (4-31) 



i=0 Ai jTOi+ i 



n 

h{U m ) = T2 \\Ri(U m )f, (4-32) 

and S(lu) is a stability factor, depending on the finite element eigenvalues {X^jLi 
and the frequency ui, defined by 



S(uj) = sup ; — — , (4.33) 
w^-c^+^c 2 



where Cj = a\j + (3. 

The terms l\ and I2 are given in matrix form by 



N 1 

h(U m ) = J2 t Rf MjRj, (4.34) 

N 1 

I 2 (U™) = £ T2 Rf MjRj. (4.35) 

i=0 i,m t +l 

Proof. As in Q3] we split the dual solution <& into two parts 4? = <&o + Here 
$0 satisfies 

a(u, * ) = (v, ip), Vv G V h , (4.36) 

and $1 satisfies 

a(v, *i) + *i) - w 2 (t>, *i) = w 2 (u, * ) - uj(Vv, * ), G V^ 1 . 

(4.37) 

Introducing this dual split in (|4.15p , using the Galerkin orthogonality prop- 
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erty (|4.1j) to subtract we obtain 

H(E) = a(E, K$o - V$ ) + tu>("DE, K<f> a - V$o) (4.38) 

-uj 2 (E,K<f>o-T<f>o) (4.39) 

+ a(E,K$i -V&i) + iU}(VE,K$i-V<f>i) (4.40) 

-lo 2 (E,11$i (4.41) 



= ^(^(C/ m ),^* - W*o) (4.42) 

i=0 

n 

+ ^(i?.f((7 m ),^* 1 -^ l * 1 ) (4.43) 

i=0 

n 

< E ||i^(L/ m )||||fti* ~ ViHi* \\ (4.44) 

n 

+ E II^^JIIHTZi*! - ViHi&xl (4.45) 

i=0 

We choose i]> = CE in the dual problem (|4.14[) . which gives (E, if)) — |||-E||| 2 . 
Using property (|4.13[) in Lemma 1 together with the Cauchy-Schwarz inequality 
on the sums in (I4.38[) . we then have 



E t — \\ R i( um )\\ 2 E m^*oiii 2 ( 4 - 46 ) 

\i=0 ! ' m,+1 / Vi=0 / 

\j=0 *>"H- 



|£||| 2 < 



\ 1/2 / „ \ 1/2 

jl^(c/-)ll 2 J (Ell^*ill 2 J • 

Using that the decomposition of V h — ©? =0 ^ is a-orthogonal, we have 

n 

Elll^*o||| 2 = |||*o||| 2 , (4.47) 

i=0 

from Parseval's identity. It furthermore holds that 

n 

EllA^* 1 || 2 <2||£* 1 || 2 , (4.48) 

2 = 

cf. e.g. Ql . 

Combining the above results, we arrive at 

/ " i \ 1/2 

lll^lll 2 ^ Er ll^(^ m )H 2 lll*o||| (4-49) 

/ " i \ 1/2 

+4Ep — iifl?(^ m )ir n^iii- (4-50) 

Vj = i,mi+l / 

We now turn to the question of stability of the dual solutions &k, k = 0,1. 
Beginning with #o in (|4.36l) we have 

a(«,* ) = (4.51) 
= a(«,£7), Vt)€^, (4.52) 
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and hence |||3»o||| = lll-E'lll- Next, the solution <I>i to (I4.37[) is given by the 
Fourier expansion 

-f-f X'j - LO 2 + LUJCj 3 ' V ' 

j = l 3 J 

where Zj, j = 1, . . . , N, is a basis of elastic eigenmodes in V h , and Cj = a\j+(3. 
Using that (C a Z, v) = X a {Z, v), for v € U'\ a > 0, we then have 



(Aj - UJ 2 + LU}Cj)(\'^ -UJ 2 ~ bLOC-j) 

(^ + ^ 2 c 2 )\(Z J ,C^o)\ 2 



? (A /l - UJ 2 ) 2 + LJ 2 C 2 
] y 3 ' 3 



(4.55) 



(w 4 +w 2 c 2 )A^ 

j (\j - W 2 ) 2 + W 2 CJ 



<^ ^\ , r, ,2^ infill 3 . ( 4 - 56 ) 



and this completes the proof. □ 

Remark 3. In Theorem[TJ we may estimate further using Young's inequality to 
obtain the subspace indicators 

Va,i = -^\\RHU m )\\ 2 + ^P^\\Rl(U m )\\ 2 , i = 0,..-,n, (4.57) 

which reflect to what degree modal truncation in each subspace influence the 
energy norm of the error in the reduced solution. 

Using these indicators we may design adaptive algorithms that automati- 
cally determines suitable refinement levels in the individual subspaces. Such 
algorithms are outlined in Section [5J along with the different numerical exam- 
ples. 

Remark 4. The matrix form of equation (|4.7p for the subspace residual R^(U m ) 
reads 

Vf MVjRj = Vf (b - KU m - tDU m + w 2 MU m ). (4.58) 
In the case of the modal basis, VfMVj = I, and 

||^([/™)|| 2 =RfR 4 (4.59) 
= (b - KtJ m - iDU m + w 2 MU TO ) T Vi (4.60) 

x Vf (b - KU m - tDU m + w 2 MU m ) 
= (b - KU m - lDTJ™ + w 2 MU m ) T W, (4.61) 
x Wf (b - KU m - tDU m + w 2 MU m ). 

The subspace residuals may similarly to the Ritz projections of the dual solution 
be approximately computed. For instance using the subspaces V h ' di \ V h ' mi , 
i = 0, . . . , n, or mass lumping techniques. 

Remark 5. In practical application we cannot evaluate S(uj) exactly as it de- 
pends on the finite element eigenvalues which in general are unknown and an 
approximation of S(uj) must be used. Such an approximation may be obtained 
using a reduced model of the eigenvalue problem. 
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Algorithm 1 



1: Start with a guess of the subspace dimensions m. 
2: Solve the problem (|3.8|) for the displacements U m . 
3: Solve the dual problem (14.14[) for <&. 

4: For each subspace V i ,m * , compute the error indicators rjj^ defined in (|4. 22[) . 
and use them together with a refinement strategy, see Remark [6l to decide 
which subspaces are eligible for refinement and how much those subspaces 
should be refined. Refine those subspaces accordingly. 

5: Repeat steps 2-4 until satisfactory results have been obtained. 



5 Numerical Examples 

In this section we apply the above developed theory in three numerical examples. 
In the first example we describe and apply an adaptive algorithm based on the 
estimate (|4.18[) in a single load case; in the second example we describe and 
apply an algorithm based on (|4.30[) in a single load case: and in the third 
example we describe and apply an adaptive algorithm based on (14.30)) for the 
computation of a series of responses when uj varies over a given range. 

In all three examples we consider the domain Q seen in Figure [TJ The 
domain is partitioned into subdomains f2i, i — 1, ... ,6, interfacing at T. Wc 
assume that the boundary is clamped at x = and stress free elsewhere. The 
reference finite element model is piecewise linear, defined on a triangular mesh 
containing approximately 7000 elements. The material constants are E = p = 1 
and v = 0.29, and the parameters a and /3 in the Raylcigh damping are chosen 
as a = /3 = 0.025. 

Example 1. We consider the single load case (w, /,grjv) where uj = 1, / = 0, 
and gjy = [0, — exp(— 100|a; — a; | 2 )], with x = (0.7,0.5). We assume that 
the goal of the computation is to control the absolute error in the functional 
H(-) — (-,TTip), where ip = [exp(— 100|a; — cci| 2 ),0], and ir is the nodal inter- 
polation operator on V . This roughly amounts to accurately computing the 
displacements in the x direction near the point Xi = (0.95,0.25). We use an 
adaptive algorithm of the form outlined in Algorithm [T] 

Remark 6. We use the following adaptive strategy: compute the normalized 
subspace indicators rji by r)i = r\ij ^2 r\i . With the objective of adding a maxi- 
mum of NMODES eigenmodes in a maximum of NITS iterations, add 

l t = [fji x NMODES/NITSJ (5.1) 

modes in subspace i = 0, . . . , n, each iteration. 

Choosing the parameters NMODES = 200 and NITS = 10 in the adaptive 
strategy (|5.1[) . the objective of the computation may be viewed as computing 
the output in H(-) as accurately as possible using approximately 200 DOF dis- 
tributed over the course of 10 iterations. Using such an objective is motivated 
by considerations regarding available precomputed basis functions and compu- 
tational resources. 

We start the algorithm with the subspace dimensions rrii = 1. Each iteration 
we compute an approximate dual solution $ using rrii + 10 eigenmodes in each 
dual subspace basis. 
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Figure 1: The domain SI partitioned into subdomains Sl^ interfacing at T. 



In TableJT] we see the subspace dimensions evolving as the adaptive algorithm 
proceeds, and in Figures we see the corresponding absolute error \H(E)\ 
together with the estimated error as the adaptation proceeds. 

For comparison we have also run the algorithm using an exact dual solution 
The resulting subspace dimensions are displayed in Tabled and the error and 
estimate can be seen in Figure [3J We see that the adaptation is similar in both 
cases, and that the estimate with approximate dual is accurate, although a slight 
underestimate is introduced after the fourth iteration. The underestimation may 
be alleviated by refining the dual more aggressively during adaptation. 

10"' 
10" 2 
10" 3 
10" 4 
10" 5 
10" 6 
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Figure 2: Semilog plot of the absolute functional error and the estimate com- 
puted using an approximate dual solution vs. the number of DOFs as the 
adaptive algorithm proceeds in Example 1. Legend: square, error; circle, esti- 
mate. 
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DOFs 

Figure 3: Semilog plot of the absolute functional error and the estimate com- 
puted using the exact dual solution vs. the number of DOFs, as the adaptive 
algorithm proceeds in Example 1. Legend: square, error; circle, estimate. 



iter. 


m 


777,1 


m 2 


777,3 


777-4 


m 5 


m 6 


2 


20 


1 


1 


1 


1 


1 


1 


4 


20 


5 


13 


10 


2 


2 


9 


6 


26 


19 


14 


14 


5 


8 


9 


8 


37 


29 


18 


15 


9 


10 


10 


10 


41 


36 


23 


15 


12 


12 


22 



Table 1: Iteration number and subspace dimensions as the adaptation proceeds 
in Example 1 using an approximate dual solution. 



iter. 


m 


777,1 


m 2 


777 3 


777,4 


m 5 


m e 


2 


20 


1 


1 


1 


1 


1 


1 


4 


20 


6 


13 


10 


2 


2 


9 


6 


26 


20 


15 


14 


3 


9 


9 


8 


33 


30 


20 


16 


8 


9 


14 


10 


35 


39 


24 


24 


8 


10 


23 



Table 2: Iteration number and subspace dimensions as the adaptation proceeds 
in Example 1 using an exact dual solution. 

Example 2. Next, we consider the load case (w, /, <7at), with lj = \/3/2, 
/ = 0, and <7at = [exp(— 100|:c — xq| 2 ), 0], with x = (0.9,0.25). In this example 
we aim to control the energy norm of the error \\\E\\\ as efficiently we can. We 
use an adaptive algorithm of the form outlined in Algorithm[2] Again we use the 
adaptive strategy ([BTTj) with the parameters NMODES = 200 and NITS = 10, 
and we start the algorithm with subspace dimensions rrii = 1. The stability 
factor S(co) is approximated using the eigenvalues from the reduced eigenvalue 
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Algorithm 2 



1: Start with a guess of the subspace dimensions m. 
2: Solve the problem (|3.8|) for the displacements U m . 

3: For each subspace V i ' m< , compute the error indicators r\ a ,i defined in (|4.57[) . 
and use them together with a refinement strategy, see Remark [6j to decide 
which subspaces are eligible for refinement and how much those subspaces 
should be refined. Refine those subspaces accordingly. 

4: Repeat steps 2-3 until satisfactory results have been obtained. 



problem: find (X m ,Z m ), such that 

a(Z m , v) = \ m {Z m , v), Wv e V h ' m . (5.2) 

We remark that although the size of the stability factor is not crucial for the 
guiding of the adaptive algorithm in this example, having a reasonable estimate 
is however important for quantitative error estimation. 

In Table [3] we see the subspace dimensions as the adaptation proceeds, and 
the corresponding error and estimate can be seen in FigureHJ From the table we 
see that the subspaces V i ' mi are refined symmetrically as should be expected 
from the given load. We see in the figure that the estimate provides an accurate 
bound on the error. 




10 ' 1 i ' ' • 1 1 ' 1 

20 40 60 80 100 120 140 160 180 

DOFs 

Figure 4: Semilog plot of the energy norm of the error and the estimate com- 
puted using the exact stability factor vs. the number of DOFs, as the adaptive 
algorithm proceeds in Example 2. Legend: square, error; circle, estimate. 



Example 3. In this example we're interested in computing the frequency re- 
sponse for a set {(w^, fk,gN,k)} of load cases. We let the goal of the computation 
be to control the relative error measured in energy norm |||_E|||/|||{7||| for each 
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iter. 


m 


mi 


m 2 


m 3 


1114 


m 5 


m 6 


2 


20 


1 


1 


1 


1 


1 


1 


4 


20 


5 


8 


8 


1 


1 


19 


6 


20 


11 


9 


9 


1 


1 


44 


8 


20 


15 


11 


11 


2 


2 


68 


10 


20 


18 


13 


13 


3 


3 


93 



Table 3: Iteration number and subspace dimensions as the adaptation proceeds 
in Example 2. 



load case, and we assume that the following estimate 

g<p^(^+SM^), (5-3) 

holds approximately. 

An adaptive algorithm designed to handle this setting is outlined in Algo- 
rithm [3] The algorithm utilizes that if a basis has been adaptively constructed 
for a given u 6 M, that basis is likely well suited for the case u + e as well, 
when e is small and that the load pattern has not changed significantly. The 
algorithm refines the subspaces V^'™ 1 ' contributing to the error and coarsens the 
subspaces that do not in order to keep the dimension of the reduced subspace 
as small as it can. 



Algorithm 3 

l: For each load case (uik, fki9N,k) in a given set. 
2: Start with a guess of the subspace dimensions rrifc. 
3: Solve the problem (I3.8[) for the displacements U™ k . 

4: For each subspace V; ' m *, compute the error indicators rj a ,i defined in (|4.57[) . 
and use them together with a refinement strategy, see Remark [71 to decide 
which subspaces are eligible for refinement/coarsening and how much those 
subspaces should be refined/coarsened. Refine/coarsen those subspaces ac- 
cordingly. 

5: Repeat steps 3-4 until satisfactory results have been obtained. 
6: Let the resulting subspace dimensions be the starting guess for the next load 
case. 



Remark 7. In Algorithm 3 we use a refinement strategy based on the following 
reasoning: begin by choosing a tolerance TOL, and let the objective for each 
load case be to refine the model such that the estimated relative error is the 
same as this tolerance, that is T] a ,i/\\\U m \\\ ~ TOL. Squaring both sides, 
we obtain 



TOL 2 . (5.4) 



\\\U m \\\ 

Assuming that each subspace should contribute equally to the error, so that 

Va.i ~ -y^Va.i, (5.5) 



14 



we have that each indicator rj a i should fulfill 



By studying the difference 



~ T0L2 f* to 

- w — — . (5.6) 



Va, TOL 2 

, z = 0,...,n, (5.7) 



we obtain a subspace indicator T 0j j, that is positive if refinement is required and 
negative if coarsening is required. By normalizing each indicator, we obtain a 
rough measure of how much each subspace should be refined or coarsened for 
(153)1 to hold true. 

Hence, let C = l/J2\ T a,i\' an d choose A i: i? 4 e N < Mj e N, where Mj is the 
number of precomputed eigenmodes for the ith subspace. Then, if r 0) i > 0, add 
LCr Qii AiJ consecutive eigenmodes are to the ith basis subject to dimlA ' mi < 
Mi, and if T a .i < 0, remove the last \Cr a ,iRi\ modes from the zth basis, subject 

tO IA > 1. 

Further, if for some load case fk>9N,k) an d subspace V^'™ 3 it holds that 

dimV/ mj = Mj, and < T a ,j < t qj , i ^ j, and v/E^/lll^lll > TOL, we 
consider the load case non resolvable given the current tolerance and maximum 
number of modes Mj, and we let the algorithm continues to the next load case. 

We let the load cases in the example be defined by u\ — 0.1, 0.2, . . . , 10, with 
/, and <7jv constant, chosen as in Example 1. The parameters in the adaptive 
strategy outlined in Remark [7] are chosen as TOL = 0.1, Ai = Mi/10, and 
Ri = Mi/10, where the number Mj of precomputed eigenmodes are Mq = 116, 
Mj = 220, i — 1,...,7. As in Example 2, we solve the reduced eigenvalue 
problem (|5.2D for a sufficiently large set of eigenvalues needed to approximate 
the stability factor 5(w). 

We start the algorithm with the subspace dimensions m< = 1. The algorithm 
terminates after requiring a total of 59 refinement iterations in order to satisfy 
the error tolerance in each of the 30 load cases. 

In Figure[5]we have plotted the energy norm of the solutions, relative errors, 
estimated relative errors, and stability factors, for each of the computed load 
cases. We see that overall the error is estimated to a high degree of accuracy, 
and we see that the estimate is close to the desired tolerance TOL=0.1, in 
every load case. The stability factor is large near ui 2 — 2.0, due to proximity 
of eigenvalues at u 2 = 1.6473 and uj 1 — 1.9812. Since the norm of the solution 
increases when approaching these values, the relative error decreases, however, 
leading to less accuracy in the estimated error. This reflects the fact that error 
estimation is more difficult near resonance frequencies. 

In Table 2] we have displayed the obtained subspace dimensions and total 
dimension, the number of required iterations, and the efficiency index EI = 
r)a,i/\\\E\\\, for each load case in the range 0.1 < lu 2 < 3.0. We see that 
typically only one or two iterations is required for each load case, except in a 
few cases. 



15 




Figure 5: Semilog plot of the energy norm of the solutions, the errors, estimated 
errors, and stability factors for each load case in Example 3. Legend: plus, 
energy norm; square, relative error; circle, estimated error; diamond, stability 
factor. 



6 Summary and Outlook 

We have presented an a posteriori error analysis for reduced finite element mod- 
els of the frequency response problem in linear elasticity constructed using com- 
ponent mode synthesis. We have derived estimates for the error in the displace- 
ments measured in a linear goal functional, as well as for the error measured 
in the energy norm. The estimate reflects to what degree each CMS subspace 
influence the error in the reduced solution allowing the design of adaptive algo- 
rithms that automatically determines suitable subspace dimensions. We have 
demonstrated our results in several numerical examples. The numerical results 
follow the theoretical predictions to a high degree of accuracy. The future of 
this research concerns its application in real world three dimensional examples. 
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Table 4: Load case, resulting subspace dimensions, total number of degrees of 
freedom, number of required iterations, and the efficiency index, for each load 
case in Example 3. 
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